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Abstract 

The quantum mechanics of two-electron systems is reviewed, starting with the ground state of 
the helium atom and helium-like ions, with central charge Z > 2. For Z = 1, demonstrating the 
stability of the negative hydrogen ion, H~ , cannot be achieved using a mere product of individual 
electron wave functions, and requires instead an explicit account for the anticorrelation among the 
two electrons. The wave function proposed by Chandrasekhar is revisited, where the permutation 
symmetry is first broken and then restored by a counter-term. More delicate problems can be 
studied using the same strategy: the stability of hydrogen-like ions (M + ,m~ ,m~) for any value 
of the proton-to-electron mass ratio M/m; the energy of the lowest spin-triplet state of helium 
and helium-like ions; the stability of the doubly-excited hydrogen ion with unnatural parity. The 
positronium molecule (e + , e + , e~, e~), which has been predicted years ago and discovered recently, 
can also be shown to be stable against spontaneous dissociation, though the calculation is a little 
more involved. Emphasis is put on symmetry breaking which can either spoil or improve the 
stability of systems. 
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I. INTRODUCTION 



The chapter on two-electron atoms or ions is of great importance when teaching quantum 
mechanics, and usually the opportunity of a transition from simple binary systems to more 
complicated structures, with examples of application of perturbation theory and variational 
methods. 

Historically, understanding the two-electrons atoms was crucial to demonstrate that the 
theory of quanta was not just an ansatz that works fortuitously for the case of the hydrogen 
atom. Indeed, while the Bohr-Sommerfeld quantization method accounts efficiently for the 
one-electron atoms, it first faced serious difficulties for the description of helium. Then 
HeisenbergP and other pioneers (for refs., see, e.g., the book by Bethe and SalpeteiP) showed 
that the helium atom can be well described in the framework of the new quantum mechanics. 

However, binding two electrons to a helium nucleus or a heavier nucleus with charge 
Z > 2 is rather obvious, as the first attached electron leaves a positively-charged kernel 
that easily traps the second electron. The problem here is to calculate approximately the 
energy spectrum and the associated wave functions and not to demonstrate the existence of 
bound states. It is thus unfortunate that many textbooks, even among the best ones, are 
restricted to the case of helium and do not discuss the more challenging case of Z = 1, i.e., 
the negative hydrogen ion. Noticeable exceptions are Refs. I3f4"f5ll6l 

The negative hydrogen ion, H~, enters a variety of physical, chemical, biological and 
geological processes.^ In astrophysics, it plays a role at Sun's surface, and its absorption 
and emission properties have been studied by Chandrasekhar in a series of papersP In some 
laboratories, there are nowadays beams of H~ . Intense beams of H~ are foreseen for future 
nuclear-fusion devicesP When teaching few-body quantum mechanics, it is the simplest 
prototype of fragile structure, at the edge between binding and non-binding, which cannot 
be described by simple tools such as Hartree wave-functions, however efficient are these 
methods to account for the properties of well-bound systems. Other examples are atomic 
clusters made of noble-gas atoms, or Borromean nuclei with two weakly-bound peripheral 
neutrons. 

Intimately related to H~ is the positronium molecule (e + , e + , e~, e~) predicted in 1945 
by Wheeler^ and discovered only recently.^ Demonstrating its stability against dissociation 
into two positronium atoms can be done with a generalization of the Chandrasekhar wave 
function, though the calculation becomes slightly more intricate. This molecule has many 
symmetries. It can be seen that breaking particle identity and breaking charge conjugation 
have dramatically different effects on its stability. In the former case, it quickly disappears, 
while it the later case, it is reinforced. In particular, the stability of the hydrogen molecule 
can be - somewhat paradoxically - demonstrated as a consequence of the stability of the 
positronium molecule. This is of course at variance with the more physical starting point of 
two infinitely massive protons, but illustrates the importance of symmetry breaking which 
enters many other fields of physics. 

This paper is aimed at reviewing what can be taught on the quantum mechanics of the 
two-electron atoms and molecules at the elementary or more advanced level. We begin in 
Sec. [IT] with the ground state of two-electron atoms and ions, which is a spin-singlet con- 
figuration. The easiest case of a central charge Z > 2 is briefly reviewed, before discussing 
the case of the hydrogen ion with Z — 1. We focus on the beautiful solution proposed 
by Chandrasekhar which is a mere product of single-electron wave functions with dif- 
ferent range parameters, supplemented by a counter-term in which the two electrons are 
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interchanged, so that the overall permutation symmetry is restored. Two other levels of 



helium-like systems are presented in Sec. Ill, the lowest spin-triplet state, whose orbital 



wave function is antisymmetric, and the unnatural-parity state of the hydrogen ion which is 



very loosely bound below its threshold. In Sec. [TV} the case of the positronium molecule is 
presented, as well as some of its less symmetric variants. After a brief summary is in Sec. [V] 
some details about the calculation of the matrix elements are given in Appendix. 

II. THE GROUND STATE OF TWO-ELECTRON ATOMS AND IONS 

We consider first the non-relativistic Hamiltonian describing two electrons of mass m and 
charge e around a fixed charged Ze, 

n 2 t> 2 Ze 2 Ze 2 e 2 
2m 2m r 1 r 2 r 12 

with ri2 = |r*2 — t*i|. The Coulomb problem has very simple scaling properties: the energies 
are proportional to e A m/h 2 ~ 27.211 eV and the distances to h 2 /(me 2 ). We shall give all 
results in natural units which correspond to treating ([!]) as if m/h 2 = e 2 = 1. The orbital 
wave function should be antisymmetric for a spin triplet, and symmetric for a spin singlet, 
as the ground state we shall consider first. 

A. The helium atom and the heavier ions 

The case of Z > 2 is treated in most textbooks. Hence we shall give only a minimal 
review, for the sake of completeness. If the last term of is omitted, the Hamiltonian is 
exactly solvable, and for the ground state, the unperturbed energy is E = -Z 2 and the 
wave function \I> z = (Z 3 /ir) exp[— Z{r\ + r 2 )]. To first order, the energy is approximated 
and upper bounded by E + E 1 = —Z 2 + z \ r i2\^ z) ( E + E\ is the variational energy 
corresponding to the trial wave function \& ' %). 

The matrix element z\f\2 \^z) is routinely estimated by a partial-wave expansion. It 
is sufficient, as done, e.g., by Peebles,^ to evoke the Gauss theorem, which states that the 
potential created at distance r 2 by a spherical shell Sqi of radius ri is 5gi/r 2 if T\ < r 2 and 
dqi/ri if r\ > r 2 . This gives 



oo 



¥■(2) 



E 1 = 4Z° /exp(-2Zr) r 2 dr / e M=^l r « df , + j exp(-2Z/) ^ <y 

Lo r 

Then for Z = 2, one obtains an energy —2.75, to be compared to E — —2.90372 . . . from the 
most sophisticated estimates, 1 14 1 15 1 and the lowest dissociation threshold E t ^ = —2. However, 
it is easily checked that this approach requires Z > 5/4 to bind two electrons. 
An easy and pedagogically instructive improvement consists of replacing ^ Q by 

# a = (a 3 /vr) exp[-a(n + r 2 )] , (3) 

where a is a variational parameter, whose value measures the effective charge seen by each 
electron. The matrix elements are the same as for a — Z, and the variational energy reads 

2 

(4) 



E = min 
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a — 2Z a H 
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the minimum being reached for an effective charge a = Z — 5/16. For Z = 2, this gives 
an improved E ~ —2.8477. Still binding with this wave function is demonstrated only for 
Z > 1.067. Thus Z = 1 requires another treatment, as described in the following section. 

B. The negative hydrogen ion 

Variational wave functions that bind H~ have been written down by Bethe, Hylleraas 
and several others. See, e.g., the book by Bethe and SalpeterP For instance, a correlation 
factor (1 + /3r 12 ) or exp(7r 12 ) can be inserted into the wave function ([3]). The most elegant 
solution is perhaps that of Chandrasekhar wri i cri reads (unnormalized) 

$ = exp(— ari — br 2 ) + e exp(— bri — ar 2 ) , e— +1 , (5) 

where the permutation symmetry is explicitly broken by a 7^ b and restored by the second 
term. As compared to the standard shell-model wave function ([3]), sometimes labelled (Is) 2 , 
this wave-function is named "unrestricted" by GoddarcP^ who gives a generalization. 

The matrix elements of $ involve the same basic integrals as for the simpler wave function 
and it is straightforward to derive, for the kinetic (T) and potential (V) energy and for 
the normalization (N), 

(&\H\$) T + V 
($|$) ~ N ' 

1 8e _ 1 1 8abe 

8a 3 b 3 + {a + bf ' ~ 16a6 3 + 16a 3 6 + (a + bf ' 

Z Z 8Ze 5e a 2 + 3ab + b 2 

~8aW ~ 8a 3 b 2 ~ (a + 6) 5 + 2(a + bf + 8a 2 b 2 {a + 6) 3 ' 

where the attractive terms (proportional to Z) are supplemented by the contribution from 
l/r 12 - 

As H~ is weakly bound, the physical picture is that of one electron far away, and the 
other one near the nucleus, feeling an almost unscreened Coulomb potential. This suggests 
the following approximation: one freezes out a = Z = 1, i.e., assumes that one of the 
electrons is unperturbed, and varies b, to get a first minimum E(l,bo) — —0.5126 that 
already establishes binding! This minimum is reached for 60 — 0.279. See Fig. [TJ dashed 
curve. 

FIG. 1: One parameter mini- 
mization of the variational en- 
ergy E(a,b) of H~, obtained 
from the Chandrasekhar wave 
function: with the approxima- 
tion of a frozen a = 1 (dashed 
curve) and, without approxima- 
tion but using the virial theo- 
rem, which removes one param- 
eter (solid curve). 



E{a, b) = 
N = 
V = 
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Using any standard minimization software easily leads to the best minimum E(ai,b\) ~ 
—0.5133 for a\ ~ 1.039 (indeed, very close to the previous approximate a = 1) and b\ ~ 0.283 
(or <2i <-> b\). For comparison, the best non-relativistic energy for an infinitely massive 
proton 

Kives™3 a bout E = -0.52775. As seen m Fig. § the minimum is not extremely 
sharp, however, the stability criterion E(a, b) < 0.5 clearly requires 6<o (or a <C b). 
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FIG. 2: Contour plot of 
the variational energy E(a,b), 
given by Eq. of the Chan- 
drasekhar wave function 
The symmetric part where a < 
b is not shown. 
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The task of minimizing E(a, b) or any similar variational energy can be simplified by 
using the virial theorem, which also holds for the best variational solution, with the mild 
restriction that the set of trial functions is globally invariant under rescaling. This was 
noticed very earlyP^^ See, also, Refs. I20f21|l22l A simple derivation of the virial theorem 
consists, indeed, to impose that in a rescaling \&(t*i, . . .) — > A -3 ™/ 2 ^(Arq, . . .), where n is the 
number of internal variables and the factor keeps the normalization, the expectation value 
of the Hamiltonian remains stationary near A = 1. This obviously works for both the exact 
solution or the best variational approximation in a given set. For instance, in the case of 



the Chandrasekhar wave function, one can set a = ao(l + x), b = ao(l 



x 



For given x, 



the minimization over the overall scale ao fixes the proper balance of kinetic and potential 
energy, as required by the virial theorem. One is left with minimizing 



V 2 



ANT 



(7) 



over the single variable x, to recover the minimum at E ~ —0.5133, as shown in Fig. [T] (solid 
curve) . 

Some years ago, another trial wave function was proposed, in a paper^l where the work 
of Chandrasekhar is not cited. It is the very compact tp = exp[— ar < — 6r>], where r< = 
min(r 1 ,r 2 ) and r< = max(r!,r 2 ). The calculation of the matrix elements of tp is similar to 
these of ([5]). Optimizing a and b with or without the help of the virial theorem leads to 
a variational energy —0.506 which demonstrates the stability of H~, but gives less binding 
than the wave function ([5]), due to the unphysical discontinuity of the radial derivatives at 
ri = r 2 . 
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C. Varying the proton charge 



The method used for H can be applied to other values of the central charge Z. One can 
first return to Z > 2, where the wave function ^ gives an energy E ~ —2.8757, instead 
of E ~ —2.8477 from the factorized wave function (pj). As 2 increases, the improvement 
becomes less and less significant, i.e., the factorized wave function ^ works almost equally 
well. See Tables HI and IPT1 



TABLE I: Binding energies (in natural units) for a series of central charge Z and electron spin S. 
The experimental energy E exv is taken from the current data bases^ES and compared to the best 
non-relativistic calculation with an infinitely massive nucleus, -E^n r! 17 * 31 - ^; the simplest Hartree- 
Fock type of calculation with an effective charge, Ef ac corresponding to ^ for S = and ( 15 ) for 
S = 1, and the energy from the Chandrasekhar wave function ^ with e = +1 for 5 = and — 1 
for S = 1, with optimized range parameters a and b. 



z s 


E 




-Efac 


Eq a b 


1 





-0.5274 


-0.5277 


-0.4727 


-0.5133 1.04 0.28 


2 





-2.9034 


-2.9037 


-2.8477 


-2.8757 2.18 1.19 


2 


1 


-2.1750 


-2.1752 


-2.1666 


-2.1607 1.97 0.32 


3 





-7.2800 


-7.2799 


-7.2227 


-7.2488 3.29 2.08 


3 


1 


-5.1103 


-5.1107 


-5.1026 


-5.0718 2.93 0.60 


4 





-13.657 


-13.656 


-13.598 


-13.623 4.39 2.98 


4 


1 


-9.2988 


-9.2972 


-9.2892 


-9.2240 3.89 0.88 
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-59.195 


-59.157 


-59.098 


-59.122 8.68 6.69 


8 


1 


-38.579 


-38.545 


-38.537 


-38.233 7.73 2.00 



Alternatively, one can investigate how far one can decrease Z without breaking the sys- 
tem. It can be checked that stability remains down to Z ~ 0.949 with the wave function (|5]). 
In the literature, the most sophisticated estimate of the critical charge to bind two electrons 
is Z c ~ 0.9107, from the so-called 1/Z expansion. By scaling, the Hamiltonian is rewritten 

as 

H _ p\ + pi _ 1 _ 1 + 1 ^ 

Z 2 2m 2m r\ r2 Z >i2 ' 

and the energy is expanded in powers of 1/Z. From g, the first terms are Z 2 [-l + 5/(8Z) + 
• ■ • ]. Elaborate studies ^ 1 25 1 26 1 (beyond the scope of the present article) have shown that 

• The expansion looks like 

5 0.157666429 0.008699032 0.000888707 \ 

E = - z { 1 -sz + — — + z> + - ■ (9) 



• The convergence of the series is associated to well isolated bound states lying below 
the threshold. When the convergence breaks, the binding is lost. The radius of 
convergence 1/Z ~ 1.098 leads to the critical charge Z c ~ 0.9107. 
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D. Proton recoil 



Coming back to unit charges, one might examine the effect of a finite mass for the proton, 



and more generally, study the stability of (M ± . 
The Hamiltonian now reads 



m T ,m~ 



") as a function of the mass ratio MJ m. 



H 



vl 



i 1 1 

2m ' 2m ' 2M r\ r 2 r\i 



(10) 



The negative hydrogen ion corresponds to M/m — > oo, the positronium ion, Ps _ , to M/m = 
1, and the molecular- hydrogen ion, H2 + , to M/m — * 0. 

Hilp3 has shown that Chandrasekhar's wave function ^ demonstrates binding for any 
M/m. When M/m is varied, the minimum is reached with the same h/a and the same 
quality of binding, as measured by the ratio of the best variational energy to the threshold 
energy 

The proof is rather straightforward. One introduces standard coordinates 



x 



ri - r 3 



y = r 2 -r 3 



R 



mr 1 + mr 2 + Mr 3 
2m + M 



(11) 



and the conjugate momenta 

(m + M)p 1 - m(p 3 + p 2 ) 

Px 



Pu 



(m + M)p 2 - m(p 3 + p^ 



M + 2m 

in terms of which the Hamiltonian becomes 



H 



4m + 2M 



+ 



Pi 
2fi 



+ 



\x\ 



2/i 



M + 2m 



1 1 

+ i 

\y\ \ x ~ 



Pl + P2 + P3 , 

(12) 



y\ 



M 



PxPy 



(13) 



Once the center-of-mass motion removed, one is left with the Hughes-Eckart term, and, in 
the bracket, a rescaled version of H~ with an infinitely massive proton and two electrons 
whose mass is decreased from m to /i — mM/(m + M). This Hamiltonian in the bracket, 
if alone, gives stability with respect to the decay into a (M + ,m~) atom and an isolated 
negative charge. However, the Hughes-Eckart term has zero expectation value within the 
wave function ^ with r\ — > \x\ and r 2 — > \y\, or in any similar wave function in which 
there is no dependence upon the angle between x and y. Hence the Chandrasekhar wave 
function gives the same energy as for the original H~, apart from an overall scaling factor 
/i/m. 



E. Symmetry breaking in three-charge systems 

So far, we studied configurations of the type (M + ,m~ ,m~) with two identical negative 
particles. One might address the question of stability of more general mass configurations 
(M + , m^, m 2 )- The most general case is discussed in the literature,^ with stable configura- 
tions such as Ps~ or H - , and unstable ones such as (p,p, e~). We shall restrict the discussion 
here to small differences between m\ and m%. 

It is known that symmetry breaking lowers the ground state. For instance, in one- 
dimensional quantum mechanics, h = p 2 + x 2 + Xx has a ground state at e = 1 — A 2 /4 
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shifted down by the odd term. More generally, if H = Hq + XH\, with Hq even and H\ 
odd under some symmetry, then the variational principle applied to H with the even ground 
state \l/o °f Hq as trial wave function, gives for the ground state E{X) < E(0) provided that 
(^ol^il^o) = 0. 

Hence, if the (M + ,rri2 ) Hamiltonian is split intd^ 

inn ^ tyy 

i/(M + ,mr,m 2 ) = H{M + ^-^-) + 1 4 2 (p? -pf) , (14) 

where /x is the average inverse mass, the ground-state energy is shifted down by the second 
term, i.e., E(M + ,7%) < E(M + , \T , But the gain is only at second order in 
m^ 1 — m^, and meanwhile, the lowest threshold decreases at first order, with E^M^m-i) < 
E 2 (M,fi) if mi < m 2 and thus /i < m 2 . Not surprisingly, the net result is that stability 
deteriorates as the two negative charges are given different masses. 

It is an interesting exercise to adapt the wave-function (|5| and to study the domain of 
stability of (M + ,m± ^m^) as a function of ni\ and m2, in the limit M — ■> oo. 



III. FIRST EXCITATIONS OF TWO-ELECTRON ATOMS AND IONS 
A. Spin-triplet ground state 

If the wave function (|5| is used with e = —1, i.e., in its antisymmetric version, it becomes 
a trial wave function for the lowest spin-triplet state. For H~, this level is unstable. Some 
results are shown in Table [T] for Z > 2. In particular, one gets E ~ —2.16064 for the 
lowest spin-triplet state of helium, to be compared to E ~ —2.17523 from wave functions 
with many parameters. Also shown in this Table is the result obtained from the simplest 
alternative wave function that comes to mind, 

lT , / v 0is(a,ri)02s(&,r 2 ) - 2s (&,ri)0i s (a,r 2 ) 
W a i,(ri,r 2 J - — , 

fl 3/2 a 3/2 ( 15 ) 

0i s (a, r) = —= exp(-ar) , 2s (a, r) = — = (1 - ar/2) exp(-ar/2) . 

7T V87T 



For a = b, it is a standard normalized shell-model wave function, and corresponds the 
fifth column of Table [TJ the results are slightly better than these the Chandrasekhar wave 
function. If the above "(Is) (2s)" wave function is used with different range parameters a 
and b (and thus with non-orthogonal individual wave functions), corresponding to a "unre- 
stricted" Hartree-Fock wave-function in the notation of Goddardp^ it gives slightly better 



results, especially for small Z. In the case of large Z, the wave function ( 15 ), with a — > b 
becomes exact. 



B. Towards a more accurate calculation 

For both the spin-singlet and the spin-triplet cases, the Chandrasekhar wave function, 
however astute, cannot describe completely the three-body ground-states such as H~. It 
can be improved by superposing more terms of the same kind. For instance, GoddarcP^ 
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TABLE II: Binding energies of H {Z = 1) and He (Z = 2) with an infinitely massive nucleus, 
obtained from the variational wave function (16). For He, we also show the two first excitations in 
the scalar sector: He* (para) with the same spin S = as the ground state, and He*(ortho) with a 
symmetric S = 1 and thus an antisymmetric space wave function. 



N 




H 


He 


He* (para) He*(ortho) 


1 


a = b = Z, c = 


-0.375 


-2.75 




1 


a = b , c = 


-0.47266 


-2.84766 




1 


a = 6, c ^ 


-0.50790 


-2.88962 




1 


a / fe, c = 


-0.51330 


-2.87566 


-2.16064 


1 


a / 6, c / 


-0.52387 


-2.89953 


-2.16153 


2 


a / 6, c / 


-0.52496 


-2.90185 


-2.14461 -2.17512 


3 


a / b, c / 


-0.52767 


-2.90328 


-2.14538 -2.17521 


4 


a / 6, c / 


-0.52771 


-2.90347 


-2.14551 -2.17522 




"Exacf'El 


-0.52775 


-2.90372 


-2.14597 -2.17523 



considered a symmetrized combination of products of Is, 2s, . . . , 5s orbitals with different 
range parameters. For the ground-state with spin singlet, he got E ~ —0.5138, which 
is a modest improvement as compared to Chandrasekhar's result E ~ —0.5133, which 
corresponds to restricting oneself to two Is orbitals. 

The most general scalar wave function depends on three variables, which can be chosen 
as the relative distances r\ = r^i = y, r2 = = x and r\2 = Hence the wave func- 
tions without explicit r 12 dependence will never approach the exact solution with arbitrary 
accuracy. 

Starting from ([5]), a natural extension is first 

\1/ = exp(— ax — by — cz) ± {a <-» b} , (16) 

and next, a superposition of such terms. With a single term (N — 1), one gets the results 
listed in Table |TT} Also shown are the improvement brought by superposing = 2, 3 or 4 
such terms. For larger N, the numerical minimization becomes delicate, and requires clever 
tools, such as stochastic search of the parameters. Frolov,^ Korobovp^ and others have 
developed a systematics of expansions based on such exponential terms and have obtained 
extremely accurate results. 



C. Hydrogen ion with unnatural parity 

Another challenging problem deals with states of unnatural parity. Consider again the 
(p,e~,e~) system, though similar considerations could be developed in the four-body case. 
If one neglects the spins and intrinsic parities, the ground state has angular momentum and 
parity + . It is the only level of below the lowest threshold H(ls) + e~ , as shown by 
HillP 
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However, the state with quantum number 1 + , i.e., unnatural parity, cannot decay into 
H(ls) + e~, at least as long as radiative corrections and spin-dependent effects are ne- 
glected. Its lowest threshold is H(2p) + e" at E th = —0.125 in natural units. It has been 
discovered^ that the lowest state of H~ with 1 + actually lies below this threshold, and the 
other calculation d 35 * 36 * 3 -^ have confirmed an energy E ~ —0.1253. The question is to find the 
most economical way of demonstrating this binding. 

The simplest wave function bearing the right quantum numbers for this state is (i is any 
projection of the vector product) 

\I/(a, b, c) = (x x y) i [exp(— ax — by — cz) + {a <-> b}] . (17) 

After angular integration, one is left with integrating a polynomial in x, y and z times an 
exponential, and the result can be deduced from a single generating function, as outlined 
in Appendix. It can be checked, after optimization of the range parameters a, b and c (or 
two of them if one uses the virial theorem) that this wave function just fails to bind the 
unnatural state of H~ . One needs a superposition, say 

y~] OLi ^(aj, b h Cj) , (18) 

i 

For a given set of range parameters {aj,6j,Q}, the coefficients oti, and the resulting energy 

are given by a generalized eigenvalue problem. Then the range parameters can be adjusted 

by standard techniques, if the number of terms is limited. If this number increases, special 

care is required, to avoid numerical instabilities. To simplify the minimization, one can 

extract the range parameters a,, bi and q from an arithmetic series a, a + (3, a + 2/5,. . ., 

allowing the possibility of equal values, e.g., bi = Ci. The minimization thus runs only on 

a and (3. If a < and (3 > 0, which helps introducing some anticorrelation among the 

two electrons, one should impose a* + ft* > 0, 6j + Cj > and q + > 0. In the case of 

H~ with 1 + , we demonstrated the stability with a few terms and thus confirmed the earlier 
results.™™ 

If one repeats the calculation in the case of Ps~ , one never reaches a variational energy 
below the Ps(2p) + e~ threshold. This confirms the conclusion by Mills, who found this state 
unbound.^ 



IV. TWO-ELECTRON MOLECULES 
A. The positronium molecule 

In 1945, Wheeler suggested a variety of new states containing positrons, which could be 
stable in the limit where internal annihilation (e + + e~ — > 7's) is neglected™ Among the 
predictions was the positronium molecule, Ps 2 , (e + , e + , e~, e~). In 1946, Ore, then at Yale, 
tried very hard to calculate this molecule, and concluded that it is likely to be unstable.^ 
However, the next year, Hylleraas and the very same Ore presented a beautiful proof of the 
stability,®' based on the wave function 

^ = exp(-ar i3 - br M - ar 24 - 6r 23 ) + {a <-> b} , (19) 

which is an obvious generalization of ([5]). All the matrix elements can be calculated 
analytically] 21 * 40 ! Some hints are given in Appendix. With a + b = 1 and a — b = (3, the 
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normalization, kinetic and potential energy read 

21 3f3 2 



33 33 - 22P 2 + 5/3 4 

71 = -r- 

16 16(1 -/3 2 ) 3 
19 21 - 18[3 2 + 5/? 4 

v = Y + ^a 







2\3 



2\2 



5f3 2 



21 - Qf3 2 + (3 4 
8(1 -f3 2 ) 3 
7 



(20) 



+ 



4/? 4 8/3 : 



+ 



/3 



2\4 



In 



/3 2 



and using the virial theorem, E = —v 2 / (4tn) should be minimized by varying f3. Though 
it does not include explicit dependence upon r\ 2 nor r%4, the wave function (19) suffices to 
establish binding at E ~ —0.5042 below the threshold for spontaneous dissociation into two 
positronium atoms, at E t ^ = —0.5. This energy has been much lowered by more and more 
sophisticated computations, 41 to reach about —0.51600. It was later realized that there are 
excited states, whose threshold is higher than two positronium atoms in the ground state, 
due to selection rules; for refs., see, e.g., Ref. H2J An indirect experimental evidence for the 
Ps2 molecule was reported recently,^ 62 years after its prediction! 



B. Other molecules 

Once the positronium molecule is shown to be stable, one might study what happens for 
other mass configurations. Though the hydrogen molecule (M + , M + ,m~ ,m~) is better de- 
scribed from the large M/m limit, i.e., the Born-Oppenheimer approximation, it is amazing 
that its stability can be understood from the M = m limit. It is also rather instructive to 
study whether or not symmetry breaking does improve binding. 

Indeed, a system (/i + , , \T , i.e., any rescaled version of Ps 2 , has many symmetries: 
exchange of the positive or the negative particles, and overall charge conjugation. 

Consider first a breaking of permutation symmetry, for simplicity, identically in the pos- 
itive and the negative sectors, keeping the average inverse mass constant. This corre- 
sponds to writing the Hamiltonian as^ 

tf(M+, m + , M- mr) = H^, (T, /x") + ( P j -pl+pj- P \) . (21) 

The second term decreases the energy of the molecule. However, the same effect is observed 
as for the three-body ion: the threshold decreases more substantially, benefitting from the 
property of two-body energies 

E 2 (M + ,M~) + E 2 {m + ,m~) <2E 2 {fi + ,ii-) . (22) 

Detailed numerical studies have shown that stability is lost for M/m > 2.2 (or m/M > 
2.2^21. 

If, instead, charge conjugation is broken, i.e., if the mass configuration becomes (M + , M + , 
m~,m~), the decomposition reads 

H(M + , M+, m-,m~)= H(fx + ,fx + , yT , //") + (p{ + pj - p 2 3 - pi) , (23) 

and again the four-body ground-state energy is lowered by the second term. Now, the 
threshold remains constant, at 

2E 2 (M + ,m-) = 2E 2 (fi + ,fi-) , (24) 

and thus the stability is improved. Indeed, the hydrogen molecule is bound by about 17% 
below the atom-atom threshold, whilst this is only about 3% for the positronium molecule. 
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V. SUMMARY 



We briefly reviewed how the stability of the ground state of H~, and the lowest spin-triplet 
state of helium can be reached with rather simple wave functions, whose matrix elements 
can be estimated by straightforward calculus. A more delicate-and less advert ised-problem 
is that of the stability of the unnatural parity states, which forces one to push further the 
variational expansion, in order to demonstrate binding. 

The main message is that the Hartree-Fock method, i.e., a factorized wave functions 
with suitable (anti-) symmetrization is extremely efficient for deeply-bound systems, but 
fails for demonstrating the binding of states at the edge between stability and spontaneous 
dissociation. This is also observed in nuclear physics: halo nuclei with a weakly bound 
external neutron and the Borromean nuclei with two external neutrons, require a dedicated 
treatment (A Borromean 3-body system is bound while its subsystem are unstable. For 
instance, 5 He = (a, n) and (n,n) are not bound, but 6 He = (a,n,n) is stable against any 
dissociation and only rely on f3 decay to disintegrate.). 

The strategy initiated by Hylleraas, Chandrasekhar and others consists of using a basis 
of functions where each term breaks permutation symmetry. The proper boson or fermion 
statistics is restored by superposing terms deduced by permutation. The same strategy 
guided Hylleraas and Ore when they derived the first proof of stability of the positronium 
molecule, and lies also underneath the most recent calculations of this system. For example, 
in their study of the positronium molecule,^ Varga and Suzuki used a basis of correlated 
Gaussians, and their own variant of the parameter search.® This method is more and more 
widely used in quantum chemistry and other few-body problems. It consists, if Xi,... x n 
are the internal variables, in describing the wave functions as superpositions of states such 
as 



ip = exp 



^2i<j a ij x i- x j + ' • • , (25) 



where the ellipses are meant for terms deduced by permutation, charge conjugation or any 
other relevant symmetry which can be explicitly enforced. The two-electron atoms seems 
the best introduction to the advanced ab-initio calculations. 



APPENDIX A: CALCULATION OF THE MATRIX ELEMENTS 

In this appendix, we give some hints to estimate the matrix elements. These of the 
wave functions ([3| and ^ have been given explicitly. Consider now its generalization 



exp(— ax — by — cz), or say, \a, b, c), whose symmetrized or antisymmetrized version (16) 
is used for spin-singlet and spin-triplet states, respectively. Here x = r 2 — r 3 , x = \\x\\, etc. 
The matrix elements are integrals over xyz dx dy dz, restricted by the triangular inequality, 
and can be all deduced from the generic function 

*>(«,/?, 7) = HI eM-*x-Py-l*)dxdydz= {a + w 4 + ^ + a) , (Al) 

\x—y\<z<x+y 



and its derivatives 



, j+j+i a<+'+*F 3 (a,/3,7) 
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For instance, the normalization of (JsJ) , besides a factor 8tt 2 due to trivial angular variables, 
reads 

(a, b, c\a, b, c) = G(l, 1, 1, 2a, 2b, 2c) , (A3) 
and any potential term is similar, e.g., 



(a, b, c|r 12 1 \a, b, c) = G(l, 1, 0, 2a, 2b, 2c) 



(A4) 



and this is easily extended to non diagonal terms, with 2a — > a + a 1 etc. 

Consider now the term p l of the kinetic energy. It is a linear combination of gradients 
with respect to the distances, 

Pi = H)[V z -Vj , (A5) 



this giving additional constant factors and y.z 
(a, b, c\p 2 \a, b, c) = (bb' + cc')(a, b, c\a, b, c) 



z 2 )/(2yz), namely, 



bb' + cd 



[G(3,0,0,a,b,c) - G(l, 2,0, a, b,c) - G(l, 0,2, a, b,c)] , (A6) 



where 2a = a + a', etc. 



For the wave function (17), some angular integrals should be done beforehand, and one 



is left with similar integrals over x, y and z. 

We now consider the four-body problem, with a wave function of the type 

\P = exp(-ar i3 - br M - cr 23 - dr 2 i) • 



(A7) 



If r 12 , r 13 and r 23 are chosen as the internal coordinates, for a scalar wave function and 
a scalar operator that do not depend explicitly on r 34 , one can work independently in 
he triangles (1,2,3) and (1,2,4) as done previously for the three-body systems, and after 
summing over the trivial angular variables, the integrals run over 



dr = r 13 ri 4 r 23 r 24 dri 2 dri 3 dri 4 dr 23 dr 2 4, 



(A8) 



A basic integral is 
F 4 (a, b, c, d, u) 



dr 12 dri 3 dri 4 dr 23 dr 24 
16 



exp(-ari 3 - 6r 23 - cr i4 - rfr 24 - ur 12 ) 



log 



(a — b) (a + b) (c — d) (c + d) |_ (a + c + u) (b + d + u) 



{b + c + u)(a + d + u) 



(A9) 



where the triangular inequalities are more easily accounted for by using the variables Si = 
i r u+i r 2i arid ti = r\i — r 2 i for % = 3, 4. All matrix elements are related to F and its derivatives. 
For instance, the normalization, first attractive term, internuclear and electronic repulsion 
of (A7) are 



n(a, b, c, d) 
v l2 (a,b,c, d) 



dF^a, b, c, d, u) 



dudadbdcdd 
d 4 F 4 (a, b, c, d, 0) 
dadbdcdd 



v 13 (a,b, c, d) 



9 4 F 4 (a, b, c, d, u) 



u=0 



dudbdcdd 
v 34 (a, b, c, d) = v 12 (a, c, b, d) , 



u=0 



(A10) 
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while for the kinetic energy of, e.g., the third particle, one gets 



(*\pl\V) = (a 2 + 6 2 )M*> - 2a%|(r 2 2 - r\, - r 2 3 ) / (2r 13 r 23 )|*> , (All) 

which can be expressed as a combination of derivatives of _F 4 . For a non-diagonal matrix 
elements between (A7) and an analogous function with range parameters a', . . . d', the co- 
efficients in the above expression become a a' + bb' and ab' + a' b, respectively, and the 
arguments of F 4 are taken to be (a + a')/2, . . . , (d + d')/2. 
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